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We investigate the possibility that the BEC-like phenomena recently detected 
on two-dimensional finite trapped systems consist of fragmented condensates. 
We derive and diagonalize the one-body density matrix of a two-dimensional 
isotropically trapped Bose gas at finite temperature. For the ideal gas, the 
procedure reproduces the exact harmonic-oscillator eigenf unctions and the 
Bose distribution. We use a new collocation-minimization method to study 
the interacting gas in the Hartree-Fock approximation and obtain a ground- 
state wavefunction and condensate fraction consistent with those obtained by 
other methods. The populations of the next few eigenstates increase at the 
expense of the ground state but continue to be negligible; this supports the 
conclusion that two-dimensional BEC is into a single state. 

PACS numbers: 03.75.Hh, 05.30.Jp, 05.70.Fh, 32.80. Pj 

Bose-Einstein condensation (BEC) cannot occur at finite temperature 
in a strictly two-dimensional trapped interacting gas in the thermodynamic 
limitP On the other hand, phenomena resembling BEC have been detected 
by experiment^ and Monte Carlo simulationJ^ on two-dimensional finite 
trapped systems. These could consist of "fragmented" condensates}^ in 
which atoms avalanche macroscopically not just into the ground state but 
into a finite set of low-energy states. Fragmented condensation has been 
studied in the context of spin-1 systems, double potential wells, and rotating 
systems with attractive interparticle interactions (the bibliography of Ref. |2 
contains a complete set of references); moreover, a mean-field analysis of a 
homogeneous three-dimensional Bose gas with repulsive interaction^^ rules 
out the possibility of fragmentation in such a system. To our knowledge, 
however, nobody has yet treated trapped 2D systems, and that is the topic 
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of this work. In what fohows we will derive and diagonalize the one-body 
density matrix — whose eigenvectors and eigenvalues respectively correspond 
to the one-body wavefunctions of the system and their populations — of a 
two-dimensional isotropically trapped Bose gas at finite temperature using 
the Hartree-Fock-Bogoliubov equations in the Hartree-Fock approximation. 
The results we obtain are consistent with condensation into a single state. 

The generalized criterion for Bose-Einstein condensation, first intro- 
duced by Penrose and OnsagerP relies on the properties of the one-body den- 
sity matrix (OBDM) of a system of identical bosons. This quantity, defined 
by the ensemble average n(r, r') = (^^(r)^'(r')), where ^' is a many-body 
field operator, has the following properties for the case of a trapped — and 
hence inhomogeneous — two-dimensional gas: (i) its diagonal, n(r) = n(r, r), 
provides the spatial density of the system; (ii) its trace is the total number 
of atoms in the gas: J cfTn{r) = N; {in) field operators are Hermitian, 
and can hence be expanded in the form ^'(r) = (j)k{r)hk; the annihilation 
(and creation) operators obey the standard commutation relations, while the 
functions (/'fc(r) form part of a complete orthonormal set. Use of the condi- 
tion {h^jhk) = fjSjk, where fj can be interpreted as the number of atoms in 
the jth state, yields 

n(x,x') = J;Ac/.^(x)c/>,(x'). (1) 

k 

Note that we have expressed the variables using the length and energy scales 
set by the trap: x = r{mu!/h)^/'^ is the position vector, and the density 
matrix and wavefunctions have been made dimensionless; in particular, n = 
hn/muj. It is straightforward to show that implies 

I A' n(x, x') 0, (x') = (x). (2) 

This expression can be interpreted as an eigenvalue equation for the operator 
J d?x' h{x,x'); the resulting eigenvectors are the one-body eigenstates and 
the corresponding eigenvalues are their occupation numbers. In the particu- 
lar case of the ideal gas, these populations are dictated by the Bose-Einstein 
distribution, which prescribes a condensation into the oscillator ground state 
at low but finite temperatures.^ The Penrose-Onsager criterion is a general- 
ization of this result: A system is said to have a Bose-Einstein condensate if 
its OBDM has a single eigenvalue of order N — identified with the condensate 
population — while the rest are orders of magnitude smaller; if more than one 
eigenvalue is of order A^, the system exhibits a fragmented condensate; if all 
of the eigenvalues are of order 1, the system is uncondensed. 

The density matrix (pQ), however, is not amenable to numerical diago- 
nalization as it stands, since it is a four-dimensional array with continuous 
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indices. In order to turn it into a diagonalizable matrix we must approxi- 
mate the integrals by discrete sums and average over all but two of these 
dimensions; this last condition makes the process feasible only for isotropic 
traps. In what follows we will diagonalize the averaged matrix 



1 



2-K 



n{x,x') = — V n(x,x')|^^Q. (3) 



The Hartree-Fock-Bogoliubov equations have been used with remark- 
able success in the study of Bose-Einstein condensation, and have been ap- 
plied to the study of two-dimensional systems.^ The model assumes an inter- 
particle interaction ^(xi — X2) = ■jS^'^^xi — X2), with 7 positive for repulsive 
interactions, and decomposes the field operator ^(x) in the form 

^ = (^) + ^ = <| + -(/J, (4) 

where the ensemble average $ is a real macroscopic wavefunction whose 
square is the condensate density no and the fluctuations about this average 
correspond to a thermal cloud of density n' = (ip^ijj). When we insert @ into 
the many-body grand canonical Hamiltonian and neglect the "anomalous" 
averages {tp^ip^) and (iptp), we obtain coupled expressions for the condensate 
and the noncondensate. The condensate is described by the Gross-Pitaevskii 
equation, 

V2|> + (A - x^)^ - 7(no + 2n')^ = 0, (5) 
and the thermal density in the Hartree-Fock approximatiorPl is given by 

n'(x) = log(l - exp[-/3(x2 + 2jfi - /i)]), (6) 

where we have introduced the inverse temperature /? = hujP/2 and the chem- 
ical potential fl = 2^/ Tito. Equations © and ©, together with the impo- 
sition that N = J d?x (no + n'), form a self-consistent set of equations that 
can be solved to study the complete thermodynamics of the trapped system 
at any temperature, with no adjustable parameters, given and 7. 

This identification also gives us an expression for the off-diagonal ele- 
ments of the density matrix. From the definition 

n(x,x') = (^t(x)^(x')) = l>*(x)l>(x') + (Vit(x)^(x')), (7) 

where we have once again neglected the anomalous averages, we obtairi^ 

1 °° 1 1 
n(x, x') = 4>(x)l>(x') + -t^y ^exp(--L{x^+x'^-2x- x') 

-ep [i(x2 + x'2)+7(n(x)+n(x'))-/i]), (8) 
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which results from neglecting the position-momentum commutator and using 
the free-particle result (x|e~^'^ |x') = e"*^''"''') /^^/47r/3. 

To solve the Hartree-Fock equations in the presence of a condensate 
we follow a self-consistent procedure: Initially, we assume that only the 
condensate is present and solve Eq. © with h' = 0. The wavefunction and 
chemical potential that result are fed into the Hartree-Fock expression for the 
density ©, which, when integrated, yields also a value for the noncondensate 
population A^'; one can then readjust the condensate population and solve 
the Gross-Pitaevskii equation that results. The process is iterated until the 
chemical potential and the populations stop changing. To solve the nonlinear 
eigenvalue problem Q we directly minimize the functionaP 

J[$] = J (fx [(V$)2 + + 27n'(x))$ + I (9) 

subject to the constraint J (Px^^ = 1. (5> = ^/NQ^.) We first set a grid 
of n fixed abscissas x^"^ that represent the radial coordinate x; in terms 
of the ordinates fj = <I>(Xj"'), J[^>] becomes a multivariate function, J = 
J[(/i, . . . , /n)^] = </[f], that can be minimized using standard optimization 
routines. The resulting function is known at the x*"' and can be evaluated 
everywhere else by interpolation. 

For each evaluation of J we need to be able to calculate both derivatives 
and integrals of functions that we know only at a few points. We have 
developed a method that uses as grid points the zeros of L^, the nth Laguerre 
polynomial; the calculation of J[f] is reduced to the matrix-vector product 



J[f] = 1 (2^"') [m] + {xf^f.f + 27n'(x5"V| + 



4, 

(10) 

(n) (n) 

where A = Yl^=i (^ttx^" ) /? Wj e^j is a normalization factor. The wy' e'^^ 



are Gauss-Laguerre weights, including the inverse of the weighting function,^ 

3 



evaluated at the grid points: w'f^ = (//^(x'"')) ^/x'"'. Finally, D is a differ- 



entiation matrix,-"^ whose elements are given by 

2L;.(4»')''*-' + L'„(4»>):.;»'-4»>- 

Now, the change of variables x^"' x^"' /b with any positive b will map the 

semi-infinite interval onto itself provided that we also rescale exp(xj"') — > 

exp(6xj"'), wj"' — > wj"'/^' ^iid I-* ^ ^l-^- This free parameter can be tweaked 
by hand in order to optimize the accuracy of the minimization process. The 
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Fig. 1. Eigenvalues (left panel) and eigenvectors (right panel) of the one- 
body density matrix of a two-dimensional interacting trapped Bose gas and 
the equivalent ideal gas. 



method just described yields self-consistent density profiles in mere seconds 
and is quite robust. 

The left panel of Fig. ^ shows the seven lowest eigenvalues of the angle- 
averaged OBDM corresponding to an isotropic two-dimensional interacting 
Bose gas in the Hartree-Fock approximation (gray bars) compared to those of 
the equivalent ideal gas (white bars). In both cases N = 1000 and T = 0.7Tc, 
where Tc is the ideal-gas transition temperature.^ The interaction parameter, 
7 = Svr X 0.0043, is arbitrary in 2D but corresponds to a 3D gas of rubidium 
atoms.l^In both cases we calculated the matrix on a 30 x 30 Laguerre grid. 

The ideal-gas results have been obtained by means of an exact expres- 
sion for the OBDM, 

oo l(3{ri—2) 

re(x,x') = — ^g-|csch2£/3((a:2+a;'2)cosh2£/3-2x-x')^ (-]^2) 

^ f=i 1 — e ^^1^ 

which is found by direct evaluation of and can be proved analytically 
to obey @ with populations prescribed by the Bose-Einstein distribution, 
fk = (exp(/3[efc — /2]) — 1)~^, where, in our system of units, = 2(2A; -|- 1); 
these values are represented in the figure by open circles, and the straight 
line depicts the distribution for continuous k. The good agreement between 
analytic and numerical results provides a useful check on our algorithms. 

The interacting-gas results clearly show the depletion of the condensate 
due to interactions and the complementary increase of the excited occupation 
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numbers. However, the ground-state population is still orders of magnitude 
greater than that of any other state (note that the scale on the y-axis is 
logarithmic). In this particular case, for example, the condensate number 
decreases from 463 to 363 but the population of the first excited state in- 
creases only from 8 to 11. Finally, the open diamond at the left shows the 
condensate number yielded by the self-consistent procedure. 

The right panel of the figure exhibits the corresponding eigenfunctions, 
which have been normalized to be unity at a; = 0. The solid lines represent 
the (spline-interpolated) interacting eigenfunctions; the ideal-gas eigenfunc- 
tions (dashed lines) have been calculated numerically but are indistinguish- 
able from the analytic wavefunctions (/!>n(x) = /^L„(x^)/y^. The open 
circles at the bottom show the Gauss-Laguerre points where the functions 
are known exactly and depict the condensate wavefunction that resulted 
from the self-consistent procedure. The wavefunctions decrease in value and 
are flattened at the center of the trap; on the other hand, the radius of the 
system becomes larger, and both condensate and noncondensate are affected. 

Our results lead us to the conclusion that, in finite 2D systems, and to 
this level of approximation, there is a phenomenon resembling BEC into a 
single state. This conclusion is supported further by Monte Carlo simulations 
that we have performed on squeezed 3D gases and whose results we intend 
to publish in the future. 
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